(C) REMARKS 



The Applicants would like to thank the Examiner for carefully reviewing the 
specification and claims. Please reconsider the application in view of the above 
amendments and the following remarks. 

1. Substitute Specification 

The Applicants noted a number of informalities in the specification, and submit herewith 
a Substitute Specification in which such informalities have been corrected. All of the 
informalities relate to preferences in terminology in the description, and the Applicants 
believe that such terminology preferences would be immediately recognized by a person 
of ordinary skill in the art. The Applicants believe that no new matter is added by the 

substitute specification. 

« 

2. Status of the Claims 

Claims 1-12 are pending in the application. Claims 1 and 7 are independent. Claims 1, 
3-7, and 9-12 stand rejected as anticipated by admitted prior art and an article authored by 
Min Lou et al. Claims 2 and 8 are objected to as reciting allowable subject matter and 
would be allowed of rewritten in independent form including all the limitations of the 
rejected base claim and any intervening claims. 

Claims 1 and 7 have also been amended to more properly recite what the Applicants 
regard as the invention. The method recited in claim 1 as amended is consistent with the 
example embodiment shown in the flow chart of Figure 2. Amended claim 7 recites a 
computer program that performs the method as recited in amended claim 1 . 

3. Claim Rejections - 35 U.S.C. § 102(b) 

Independent claims 1 and 7 have been amended to recite the subject matter of claims 2 
and 8, respectively. Because claims 2 and 8 were found to recite allowable subject 
matter, the rejection for anticipation is believed to be overcome. 
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The Applicants believe this Reply to be fully responsive to each and every ground of 
rejection or objection cited by the Examiner in the Office Action of June 21, 2005 and 
respectfully request early favorable action on this application. 



Respectfully submitted, 



Date: 





"Richard A. Fagin, Reg. No. 39,182 
P.O. Box 1247 
Richmond, TX 77406-1247 

Counsel of Record 

E. Eugene Thigpen, Reg. No. 27,400 

PGS Americas, Inc. 

15150 Memorial Drive 

Houston, Texas 77079 

Telephone: (281) 509-8000 
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Stable estimation of anisotropic 
parameters for p-wave prestack imaging 

Cross-reference to related applications 

Not applicable. 

Statement regarding federally sponsored research or development 

Not applicable. 

Background of Invention 

Field of the Invention 

[0001] The invention relates generally to the field of seismic data processing methods. 

More specifically, the invention relates for—to methods for estimating seismic 
compressional wave anisotropy in the velocity for use in pre-stack seismic migration 
imaging. 

Background Art 

[0002] Seismic surveying is used to evaluate structures of 5 compositions of, and fluid 

content of subsurface earth formations. A particular application for seismic surveying is 
to infer the presence of useful materials, such as petroleum, in the subsurface earth 
formations. Generally, seismic surveying includes deploying an array of seismic sensors 
at or near the earth's surface, and deploying a seismic energy source near the sensors also 
at or near the surface. The seismic energy source is actuated and seismic energy 
emanates from the source, traveling generally downwardly through the subsurface until it 
reaches one or more acoustic impedance boundaries. Seismic energy is reflected from 
the one or more impedance boundaries, where it then travels upwardly until being 
detected by one or more of the sensors. Structure and composition of the subsurface is 
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inferred from the travel time of the seismic energy, and the amplitude and phase of the 
various frequency components of the seismic energy with respect to the energy 
emanating from the seismic source. 

[0003] In order to infer the structures of subsurface earth formations from seismic travel 

times measured at the earth's surface from a source position at the surface, it is necessary 
to determine the velocity of the various formations through which the seismic energy 
passes. Velocities of the earth formations can vary both with respect to depth in the earth 
(vertically), and with respect to geographic position (laterally). Seismic data, however, 
are recorded only with respect to time. Methods known in the art for estimating 
velocities of the earth formations both vertically and laterally rely on inferences about the 
travel path geometry of the seismic energy as it travels from the seismic source to the 
various seismic receivers deployed at the earth's surface. 

[0004] Migration is a process performed on seismic data in which depth estimates to one 

or more reflective horizons (acoustic impedance boundaries) in the earth are made from 
the "two way" travel time of seismic energy from the source, to the reflective horizons, 
and back to the seismic receivers. The depth estimates are computed and displayed with 
respect to geographic position of the reflective horizons. Depth estimates based on two 
way travel time must be corrected for the effects of seismic energy travel path differences 
between various seismic energy source and receiver geographic positions that are used 
during data acquisition. In order to correct the depth estimates for source and receiver 
positions, it is necessary to accurately estimate the velocity of seismic energy in the earth 
from the earth's surface te-(or the ocean bottom in marine seismic data) to each 
subsurface reflective horizon. Methods are known in the art for estimating velocity from 
two wa[y travel time from the seismic source to the reflective horizons and back to the 
seismic receivers. 

[0005] It is known in the art that some earth formations are anisotropic with respect to 

seismic velocity. Anisotropy can be observed as different velocity in a single medium 
depending on the direction of propagation of the seismic energy through the medium. In 
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anisotropic earth formations, the depth estimating (imaging) process needs to take the 
anisotropy into account in order to accurately position the seismic reflective horizons 
with respect to depth and geographic position. Reliable estimation of anisotropy 
parameters is thus important to accurate depth imaging. Because the anisotropy 
parameters and velocity are related, however, estimation of anisotropy parameters is not 
trivial and is often unstable using methods known in the art. 

[0006] It is known in the art to estimate certain anisotropy parameters (called "anelliptic" 

parameters) from prestack seismic data. "Prestack" seismic data are seismic data traces 
(traces being a term used in the art for a record of the signal detected by one of the 
seismic receivers) that have not been summed, average, or otherwise processed together 
with any other seismic data trace. See for example, Alkhalifah, T., Velocity analysis 
using nonhyperbolic moveout in transversely isotropic media, Geophysics 62, 1839-1845 
(1997). See also Lou, M., Pham, L.D. and Lee, S., Anisotropic parameter estimation 
from joint P- and C-wave data, 64 th Ann. Mtg., Eur. Assn. Geosci. Eng., Florence (2002). 
See also, Martinez, R.D. and Lee, S., 2002, A strategy for anisotropic P-wave prestack 
imaging, 72 nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 149-152 
(2002). 

[0007] A method disclosed in the Martinez and Lee reference, for example, combines 

prestack time migrations and prestack depth migrations for estimating the anisotropic 

parameters which are related to the a nd-velocities. The method disclosed in the Martinez 

and Lee reference has been shown to generate an accurate result for -to estimate the 

effective anelliptical (rj e ff ) parameter. Figure 1 shows a flow chart of the method 

disclosed in the Martinez and Lee reference. At 10, seismic data can be prestack time 

migrated. At 11, the prestack time migrated data may be scanned to obtain a model of 

NMO velocity and effective anelliptical parameters. The NMO velocity and the effective 

anelliptical p arameter models in time may be processed, at 12, to obtain interval NMO 

velocity and anelliptical parameter models in depth . Sonic logs, check shot surveys, 

and/or vertical seismic profile (VSP) data, at 13, may be used to derive the vertical 

velocity model. Alternatively, at 16, the interval NMO velocity and the anelliptical 
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parameter models may be put through an anisotropic ray tracing procedure, and at 15, the 
data are prestack depth migrated. The interval velocity (NMO and vertical) and 
anelliptical parameter models may be adjusted in depth, at 14, and the ray tracing and 
depth migration are repeated until a satisfactory imaging result is obtained. 

[0008] The effective parameters 7j eff and V nmo with respect to two-way seismic travel 

time, obtained from the prestack time migration, need to be converted into "interval" 
parameters rj in t and V int with respect to depth, as explained in the Martinez and Lee 
reference. t] int is described in the Martinez and Lee reference as being inverted from the 
effective anelliptical parameter (rjeff). The effective anelliptical parameter is equivalent to 
a n root mean square (RMS) average^of anelliptical parameter from the earth's surface or 
water bottom to a selected depth in the subsurface ) _and e ffective velocity (JW ). -Both 
the e£feem*e~ NMO velocity and effective anelliptical parameter are typically determined 
by "scjanning" compressional (P-wave) moveout on common-image gathers from 
prestack time migration. Scanning refers to the process of calculating a correlation or 
semblance between selected seismic traces after adjusting using various values of NMO 
compressional velocity and the effective anelliptical parameter. Because the selected 
traces typically represent different source to receiver distances (offsets), the semblance 
between traces will be related to the mov e out N MO velocity and anelliptical parameters 
of the earth formations between the seismic source and the one of the seismic receivers 
used to record each trace. The values of compr e ssional N MO velocity and effective 
anelliptical parameters which provide the greatest value of semblance are selected as the 
effective (or normal mov e out) N MO velocity and effective anelliptical parameters. The 
values of effective anelliptical parameter and normal moveout N MO velocity are then 
used to obtain values of interval NMO velocity and interval anelliptical parameter. The 
inversion to obtain each set of interval values is performed using the well known Dix 
equation. The parameters determined using the Dix equation are then used to generate a 
prestack depth migrated image. The image is checked for error, and the values of interval 
NMO velocity and anelliptical p arameters are adjusted, and the prestack depth migrated 
image is generated again, as explained above with respect to Figure 1. The foregoing 
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procedure is repeated until the image error is determined to be at a minimum, or below a 
selected threshold. Image error is determined by lack of "flatness" (consistency in two 
way time depth to a particular reflective horizon after prestack depth migration) in 
common image gathers (CIGs). 

[0009] Because of its differential nature, however, the Dix equation boosts errors present 

in the inputs (the NMO effectiv e velocities and effective anelliptical parameters), making 
its output unstable. This is especially true in the inversion for the anelliptical parameter 
Tii nt because the NMO velocity is related to the anelliptical parameter by the 4 th power of 
the velocity. Thus, relatively small errors in the NMO velocity can be translated into 
much larger errors in the output interval anelliptical parameters 7] int . There thus exists a 
need for a better technique to determine interval NMO velocities and the anelliptical 
parameters from seismic data. 

Summary of Invention 

[0010] One aspect of the invention is a method for determining interval velocities and 

anelliptical parameters. The method includes determining normal mov e out NMO 
velocities and effective anelliptical parameters from seismic data traces. The normal 
moveout N MO velocities are processed to obtain interval velocities. The normal 
moveout velocities, effective anelliptical parameters and interval v e locitie s are inverted 
to obtain the interval anelliptical parameters. In one embodiment, the inversion is 
performed using preconditioned damped least squares. In one embodiment, the normal 
moveout N MO velocities and effective anelliptical parameters are determined by scanning 
a time migrated trace gather that accounts for ray bending. 

[0011] Another aspect of the invention is a computer program stormed [stored] in a 

computer readable medium. A program according to this aspect of the invention includes 
logic operable to cause a programmable computer to perform the following process. 
First, normal mov e out N MO v elocities and effective anelliptical parameters are 
determined from seismic data traces. The normal mov e out NMO velocities are processed 
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to obtain interval velocities. The normal moveout velociti e s, effective anelliptical 
parameters , NMO velocities and interval NMO velocities are then inverted to obtain the 
interval anelliptical parameters. In one embodiment, the inversion is performed using 
preconditioned damped least squares. In one embodiment, the normal mov e out N MO 
velocities and effective anelliptical parameters are determined by scanning a time 
migrated trace gather that accounts for ray bending. 

[0012] Other aspects and advantages of the invention will be apparent from the following 

description and the appended claims. 

Brief Description of Drawings 

[0013] Figure 1 shows a flow chart of prior art processing to obtain an e lliptical 

anisotropic p arameters. 

[0014] Figure 2 shows a flow chart of one embodiment of a process according to the 

invention. 

[0015] Figure 3 shows a comparison of interval anelliptical parameters obtained using a 

method according to the invention with parameters obtained using the Dix equation. 

[0016] Figures 4a and 4b show, respectively, image gathers made using a method 

according to the invention and by- using the Dix equation. 

[0017] Figures 5 a and 5b show, respectively, input values of NMO velocity and 

calculated values of interval velocity obtained using a method according to the invention 
on a seismic data survey line. 

[0018] Figure 6a and 6b show, respectively, input effective anelliptical parameters and 

calculated interval anelliptical parameters obtained for the same survey line as shown in 
Figure 5. 

[0019] Figure 7a shows part of a final image made from anisotropic prestack depth 

migration using the interval NMO velocit y, vertical velocity constrained by the well 
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depths and interval anelliptical parameters obtained using a method according to the 
invention. 

[0020] Figure 7b shows the corresponding image made from isotropic prestack depth 

migration only using the interval NMO velocity obtained using the Dix equation. 

[0021] Figure 8 shows a programmable computer used to read a program stored in a 

computer readable medium. 

Detailed Description 

[0022] Methods according to the invention are applicable to transverse isotropy with 

vertical axis of symmetry (VTI) layered earth formations. The formations also have 
vertical velocity which changes with depth beneath the earth's surface. A method 
according to one embodiment of the invention includes using a preconditioned damped 
least squares (DLS) inversion process to determine interval NMO velocities, V inh and 
interval anelliptical parameters, ^r^from the effective values of the anelliptical 
parameter and from effective normal mov e out NMO velocities. DLS inversion has been 
successfully used, for example, to invert interval velocity from RMS velocity. See, for 
example, Claerbout, J.F. and Fomel, S., 2001, Geophysical Estimation by Examples: 
Environmental soundings image enhancement, Stanford University (2001). See also 
Clapp, R.G., Sava, P. and Claerbout, J.F., Interval velocity estimation with a null-space, 
SEP 97, 147-156(1998). 

[0023] In a method according to the invention the Vi nt and 7j int inversions can be 

formulated as a system of linear equations. Therefore, convergence to a solution set of 
interval NMO velocities and anelliptical parameters (a "solution model") can be 
performed by means of conjugate gradients. In one embodiment, preconditioning is used 
to speed up convergence to the solution model. In the description which follows, the 
mathematical basis for determining interval anelliptical parameters and interval NMO 
velocities will first be provided, and then a description of one implementation of a 
method according to the invention will be provided. 
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[0024] 1. The relationship between effective and interval values of the anelliptical 

parameter 7] 

[0025] Referring, for example, to Alkhalifah, T., Velocity analysis using nonhyperbolic 

moveout in transversely isotropic media, Geophysics 62, 1839-1845 (1997), the effective 
anelliptical parameter, rj e ff 9 can be defined by the expression: 

[0026] M , 0 ) ,ILJ_ j; + ZrjMVr -l) (1) 

[0027] The above expression is related to the normal moveout equation for 

compressional (P) waves in layered VTI media, in which t 0 represents the two way travel 
time at zero offset (a condition where the seismic source and receiver are collocated at 
the earth's surface or water surface), and ^represents the interval NMO velocity, 

[0028] The normal mov e out NMO velocity (V nm0 ) can be determined by the expression: 

[0029] OO-fftfM* < 2 > 

[0030] The effective anelliptical parameter,/^, can be split into two components: 

[0031] 7* =7,+*. (3) 

[0032] In equation (3), the component r\ v can be determined by the expression: 

[0033] rj v =\\-±r f YiMrfr-ll (4) 

V, 0 nmo J 

[0034] and the component r\ a can be determined by the expression: 

[0035] tj a = -L- £° VlMHMdr (5) 

*(T nmo 



[0036] The-^ represents the effective component of rj e ff is-due to ray bending. Ray 

bending results from vertical seismic velocity variation with depth. The 7j a component 
represents the effective intrinsic anisotropy of the earth formations. rj e ff can be 
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determined by scanning gathers generated in prestack time migration that account for the 
ray bending. Therefore, t] e ff determined by scanning such gathers is essentially equal to 
7] a , the effecti ve intrinsic anelliptical parameter. This means that an inversion for rj int can 
be performed based on equation (5) using the // e # parameters determined by scanning a 
prestack time migrated gather that accounts for ray bending. 

[0037] Preconditioned DLS inversion process to obtain values of rj int 

[0038] Figure 2 shows a flow chart for an embodiment of DLS inversion to obtain 7j int 

from the interval NMO velocities, normal mov e ont NMO velocities and the effective 
intrinsic anelliptical parameter. First, at 20, seismic data are scanned to obtain values of 
V nmo , and are scanned at 22 to obtain r\ eff . Scanning, as previously explained, includes 
calculating semblance between selected seismic traces with respect to values of time, 
normal moveout N MO velocity and effective anelliptical parameter. Values of effective 
anelliptical parameter and NMO normal moveout velocity which result in the highest 
semblance values between selected traces are selected as the values of effective 
anelliptical parameter and normal moveout N MO velocity at each selected time. 

[0039] Next, at 24, the V nm0 values are inverted using DLS to obtain a stable set of values 

of V inh at 26. The stable set of values for V int 26 is important to obtain stable estimates of 
rjinu because, as can be inferred from equations (1) and (5) above, the velocity is related 
to the anelliptical parameter by the 4 th -power of the velocity. Thus, a small error in 
IW>_may be translated into large errors in the output rj int . The inverted V int 26 can also be 
used in a time to depth conversion, at 31, to initialize as part of an interval NMO velocity 
model for depth imaging, as shown in Figure 2 at 32. The effective intrinsic anelliptical 
parameter rj a is determined, as explained above, from the effective anelliptical parameter 
Hetfby using a time migration analysis that accounts for ray bending, in which case the 
effective anelliptical parameter is equal to the effective intrinsic anelliptical parameter. 

[0040] At 28, the interval anelliptical parameter is DLS-inverted from rj ay V nm(h and V int 

to generate a set of values of rj int at 29. A time-to-depth conversion can be applied, at 30 
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to the resultant 7j int values, using V inh in order to obtain its value versus depth, which can 
then be used to initialize the rj m t model for depth imaging, as shown at 34. rj a , as 
explained above, is equal to the effective anelliptical parameter when the effective 
anelliptical parameter is scanned from time migrated image gathers that account for ray 
bending. V int can also be used to generate an initial interval depth migrat e d N MO 
velocity4mag e in depth , as shown at 32. 

[0041] Having explained a method according to the invention in general terms, an 

embodiment of DLS will be explained. The objective function fitting goal of the DLS 
inversion for either Vua or T] int can be expressed as: 

[0042] (Cm - d|| + s ||Dm|| = min (6) 

[0043] In equation (6), d is the data vector, whose z th component is (V nmo )^(t 0 ) i and 

(LJ/fe)/^)! for the v int ^d rjint inversions, respectively, m is the model vector, 
whose I th component is (V int )] (square of the i th interval NMO velocity) and (V int )*(tl int )i 
(product of the fourth power of the ith interval NMO velocity and the i th interval 
anelliptical parameter) for the V int and 7] int inversions, respectively, s is the damping 
factor. C is the causal time-integration matrix, the non-zero component of which in its / th 
column is equal to (t 0 ) M -(t 0 ) r D is the differentiation matrix and is equal to the 
inverse of the causal time-integration matrix: D = C~\ 

[0044] The first part on the left-hand side of equation (6) is the data fitting objective goal , 

which is a user selected error threshold between the modeled and the actual NMO 
velocity and effective anelliptical parameter fields. The second part is the model 
damping objectiv e goal , which is selected to minimize the variability of the inverted 
model. ■ The system of linear equations represented by equation (6) can be solved by the 
conjugate gradient method, which is known in the art. 

[0045] To speed up the solution convergence, in some embodiments a preconditioning 

step can be performed. See the Claerbout and Fomel reference cited above for an 
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example of preconditioning. Preconditioning can be performed by changing variables in 
equation (6) as shown: 

[0046] p = Dm (7) 

[0047] where p is a preconditioning vector. Thus, from equation (6) a final system of 

linear equations to be solved can be expressed as: 

[0048] ||c 2 p-d|| + 4p|| = min (8) 

[0049] A method according to the invention generates a stable result for the interval 

NMO velocity and anelliptical paramete r and th e interval velocity . -The method 
according to the invention systematically takes the data fitting and model smoothing into 
account in the sense of least squares. While it would be possible to filter or smooth 
unstable interval anelliptical parameters and interval v e locity parameters v elocities 
obtained using the Dix equation, the results so obtained have been shown not to fit actual 
seismic data very well, as compared with the parameters obtained using the method 
according to the invention. The method of the invention is robust and effective, as has 
been shown by examples in this description and in actual use. 

[0050] To illustrate the performance of the method of the invention, presented here is a 

field data example. Figure 3 shows the inverted 7] inh as shown by curve 40, from the 
input values of 7] a , as shown by curve 42 at one image location. For comparison, Figure 
3 also shows the corresponding result obtained by the Dix equation, as shown by curve 
44. It is apparent from the graph in Figure 3 that the DLS approach gives a much more 
stable result. 

[0051] Figures 4a and 4b show two migrated common-image gathers (CIGs) from depth 

imaging, which were generated using the V int and rj int parameters obtained using the 
method of the invention (Figure 4a) and using the Dix equation (Figure 4b). For ray 
tracing, it proved necessary to apply smoothing to the parameters obtained with the Dix 
equation. However, it was not necessary to apply smoothing to the parameters obtained 
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with the method of the invention By using the velocity mod e l p arameter models obtained 
with the method of the invention, depth migration produced a better image. 

[0052] Figure 5 a shows a three dimensional graph of NMO velocity with respect to time 

in a selected number of seismic data traces in a survey "line." Time is on the vertical 
axis, trace position is shown in the horizontal axis, and NMO velocity is shown by the 
gray scale. The inverted interval NMO velocities obtained using the method according to 
the invention for a selected seismic data survey "line" are shown graphically in Figure 
5b, wherein interval NMO velocity is shown by the gray scale. The corresponding 
effective and interval anelliptical parameter values are shown graphically in Figures 6a 
and 6b, respectively. Anellpitical parameter values are shown by the gray scales in 
Figures 6a and 6b. Although smoothing was not applied to the parameter fields shown in 
Figures 5b and 6b, the effective and interval anelliptical p aramete r sections are shown to 
be substantially stable and smooth in both the v e rtical and lat e ral directions . 

[0053] Figure 7a shows a part of a final image, which was generated from anisotropic 

prestack depth migration using V int and rj int obtained using the method of the invention. 
Figure 7b shows the same part of the image, which was generated from isotropic prestack 
depth migration using V int obtained using the Dix equation. Compared with the image in 
Figure 7b, the image in Figure 7a is a-much better focused as illustrated in the encirciled 
areas indicated by reference numeral 50 in Figure 7a and reference numeral 52 in Figure 
7b. The image in Figure 7a also appears to correspond better to geologic information 
provided by a projected wellbore, shown at curve 46 and line position 48 in Figure 7. It 
is shown that the method of the invention can provide reliable parameters for anisotropic 
migrations. 

[0054] In another aspect, the invention relates to computer programs stored in a computer 

readable medium. Referring to Figure 8, the foregoing process as explained with respect 
to Figure 2, can be embodied in computer-readable code stored on a computer readable 
medium, such as floppy disk 68, CD-ROM 70 or magnetic hard drive 66 forming part of 
a general purpose programmable computer. The computer, as known in the art, includes 
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a central processing unit 60, a user input device such as a keyboard 62 and a user display 
64 such as a flat panel LCD display or cathode ray tube display. According to this aspect 
of the invention, the computer readable medium includes logic operable to cause the 
computer to execute steps as set forth above and explained with respect to Figure 2. 

[0055] While the invention has been described with respect to a limited number of 

embodiments, those skilled in the art, having benefit of this disclosure, will appreciate 
that other embodiments can be devised which do not depart from the scope of the 
invention as disclosed herein. Accordingly, the scope of the invention should be limited 
only by the attached claims. 
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METHOD FOR STABLE ESTIMATION OF ANISOTROPIC 
PARAMETERS FOR P-WAVE PRESTACK IMAGING 

Abstract 

A method is disclosed for determining interval anelliptical anisotropic p arameters. The 
method includes determining normal moveout (NMO) velocities and effective anelliptical 
parameters from seismic data traces. The normal moveout NMO velocities are processed to 
obtain interval NMO velocities. The normal moveout NMO velocities, effective anelliptical 
parameters and interval NMO velocities are inverted to obtain the interval anelliptical 
parameters. In one embodiment, the inversion includes damped least squares. In one 
embodiment, the inversion is preconditioned. 
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